基礎数理物理学レポート

川口健三郎

190-s1490

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os
from matplotlib.animation import ArtistAnimation
from IPython import display
from IPython.display import HTML
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
# %matplotlib inline

# plt.rcParams['animation.ffmpeg_path'] = '/opt/homebrew/bin/ffmpeg' # It's needed to convert animation plot to html5 video

関数諸々

# Lax-Wendroff two step method
def lax_w_2step(u, calw):
    # Firts step
    w = calw(u)
    u_first = np.empty_like(u)
    u_first[:,0:-1] = u[:,0:-1] - kap * (w[:, 1:] - w[:,0:-1]) 
    u_first[:, -1] = u_first[:, 0] # BdC
    # Second step
    w = calw(u_first)
    u_second = np.empty_like(u)
    u_second[:, 1:] = 0.5 * (u[:, 1:] + u_first[:, 1:] - kap * (w[:, 1:] - w[:, :-1]))
    u_second[:, 0] = u_second[:, -1] # BdC
    return u_second

def calw_fluid(u):
    w = np.empty_like(u)
    rho = u[0].copy()
    rhov = u[1].copy()
    e = u[2].copy()
    v = rhov / rho
    p = (gamma - 1) * (e - rhov * v / 2)
    w[0] = rhov
    w[1]= rhov * v + p
    w[2] = (e + p) * v
    return w


# Lax-Wendroff two step method for 2-dimention
def lax_w_2step_2d(u, calw):
    # Firts step
    w1, w2 = calw(u)
    u_first = np.ones_like(u)
    u_first[:,:-1, :-1] = u[:,:-1, :-1]\
                        - kap1 * (w1[:, :-1, 1:] - w1[:,:-1, :-1])\
                        - kap2 * (w2[:, 1:, :-1] - w2[:,:-1, :-1])
    u_first[:, :, -1] = u_first[:, :, 0].copy() # BdC
    u_first[:, -1, :] = u_first[:, 0, :].copy() # BdC
    # Second step
    w1, w2 = calw(u_first)
    u_second = np.ones_like(u)
    u_second[:, 1:, 1:] = 0.5 * (u[:, 1:, 1:] + u_first[:, 1:, 1:]\
                        - kap1 * (w1[:, 1:, 1:] - w1[:, 1:, :-1])\
                        - kap2 * (w2[:, 1:, 1:] - w2[:, :-1, 1:]))
    u_second[:, :, 0] = u_second[:, :, -1].copy() # BdC
    u_second[:, 0, :] = u_second[:, -1, :].copy() # BdC
    return u_second

def calw_fluid_q4v2(u):
    w = np.empty_like(u)
    rho = u[0].copy()
    rhov = u[1].copy()
    e = u[2].copy()
    v = rhov / rho
    p = (gamma - 1) * (e - rhov * v / 2)
    w[0] = rhov
    w[1]= rhov * v + p
    w[2] = (e + p) * v
    w[3] = w[2].copy()
    return w, w

def calw_fluid_2d(u):
    w1 = w2 = np.empty_like(u)
    rho = u[0].copy()
    rho_vx = u[1].copy()
    rho_vy = u[2].copy()
    e = u[3].copy()
    
    vx = rho_vx / rho
    vy = rho_vy / rho
    v2 = vx ** 2 + vy ** 2 # velocity squared
    p = (gamma - 1) * (e - rho * v2 / 2)
    w1[0] = rho_vx
    w1[1] = rho_vx * vx + p
    w1[2] = rho_vx * vy
    w1[3] = (e + p) * vx
    
    w2[0] = rho_vy
    w2[1] = rho_vy * vx
    w2[2] = rho_vy * vy + p
    w2[3] = (e + p) * vy
    
    return w1, w2

課程A

A1 課程1 　　拡散方程式

例題

$$ \frac{∂u}{∂t} = \frac{∂^2u}{∂x^2} \\ u(0,t) = u(1, t) = 0 \\ u(x, 0) = \sin(\pi x) $$

差分方程式は $$ u_i^{n+1} = u_i^n + \frac{s}{h}(u_{i+1}^n - 2u_i^n + u_{i-1}^m) \\ (n = 0, 1, 2, \dots, i = 1, 2, \dots, I -1) $$

# 代入とコピーでは意味合いが違うので要注意
imax = 20
x = np.linspace(0, 1, imax + 1)
y = x.copy()
z = x
x += 2
print(f"x = \n{x}")
print(f"y = \n{y}")
print(f"z = \n{z}")

x = 
[2.   2.05 2.1  2.15 2.2  2.25 2.3  2.35 2.4  2.45 2.5  2.55 2.6  2.65
 2.7  2.75 2.8  2.85 2.9  2.95 3.  ]
y = 
[0.   0.05 0.1  0.15 0.2  0.25 0.3  0.35 0.4  0.45 0.5  0.55 0.6  0.65
 0.7  0.75 0.8  0.85 0.9  0.95 1.  ]
z = 
[2.   2.05 2.1  2.15 2.2  2.25 2.3  2.35 2.4  2.45 2.5  2.55 2.6  2.65
 2.7  2.75 2.8  2.85 2.9  2.95 3.  ]

imax = 20
dx = 1/imax
x = np.linspace(0, 1, imax + 1)

# step 2
u = np.sin(x * np.pi)
plt.plot(x, u, label=f"u0")

# step 3
dt = 1e-3
kap = dt/dx**2
nmax = 70
for n in range(nmax):
    u0 = u.copy()
    u = np.array([u[i] + kap * (u[i+1] - 2*u[i] + u[i-1]) for i in range(1,imax)])
    u = np.append(0, u)
    u = np.append(u, 0)
    if (n+1) % 10 == 0:
        plt.plot(x, u, label=f"t = {dt * n:.2f}")
plt.legend()
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.plot([0.5,0.5])

[<matplotlib.lines.Line2D at 0x11f77d3d0>]

数値積分の簡単化

ちなみに、pythonではiについては次のようにfor文を用いずに計算をすることもできる。

下の関数でステップを一回進めることができる。（よって、nについてのループは必要。）

なお、bd0, bd1はx=0,1での境界条件。

# Finite-Difference Method(名前があっているかわからないけど)
def fdm(u, bd0=0, bd1=0):
    w = np.empty(len(u))
    w[1:-1] = u[1:-1] + kap * (u[2:] - 2*u[1:-1] + u[:-2])
    w[0] = bd0
    w[-1] = bd1
    return w

imax = 20
dx = 1/imax
x = np.linspace(0, 1, imax + 1)

# step 2
u = np.sin(x * np.pi)
plt.plot(x, u, label=f"u0")

# step 3
dt = 1e-3
kap = dt/dx**2
nmax = 70
for n in range(nmax):
    u = fdm(u)
    if (n+1) % 10 == 0:
        plt.plot(x, u, label=f"t = {dt * n:.2f}")
plt.legend()
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.plot([0.5,0.5])

[<matplotlib.lines.Line2D at 0x12f3a0c70>]

かなりすっきりとしたコードになり、再利用もできるため関数化は結構便利。

試してみよう2

第4章例題4.3（シルクハット型関数の初期条件）

$$ \frac{∂u}{∂t} = \frac{∂^2u}{∂x^2} \\ u(0,t) = u(1, t) = 0 \\ u(x, 0) = \Omega(x) = \begin{cases} 1 & \left( \frac{1}{3} \le x \le \frac{2}{3} \right) \\ 0 & (\text{others}) \end{cases} $$

imax = 20
dx = 1/imax
x = np.linspace(0, 1, imax + 1)

# step 2
u = np.zeros(imax + 1)
u[int(0.25/dx):int(0.75/dx + 1)] = 1
plt.plot(x, u, label=f"u0")

# step 3
dt = 1e-3
kap = dt/dx**2
nmax = 70
for n in range(nmax):
    u = fdm(u)
    if (n+1) % 10 == 0:
        plt.plot(x, u, label=f"t = {dt * n:.2f}")
plt.legend()
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.plot([0.5,0.5])

[<matplotlib.lines.Line2D at 0x12eb79dc0>]

数値不安定性

imax = 20
dx = 1/imax
x = np.linspace(0, 1, imax + 1)

# step 2
u = np.zeros(imax + 1)
u[int(0.25/dx):int(0.75/dx + 1)] = 1

# step 3
dt = 2e-3
kap = dt/dx**2
nmax = 3
for n in range(nmax):
    plt.plot(x, u, label=f"t = {dt * n:.3f}")
    u = fdm(u)
plt.plot(x, u, label=f"t = {dt * (n + 1):.3f}")
plt.legend()
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.plot([0.5,0.5])

[<matplotlib.lines.Line2D at 0x12f389070>]

魔女の帽子型関数

初期条件は $$ u(x,0) = \Lambda(x) = \begin{cases} 1 - 4|x -\frac{1}{2}| & (|x - \frac{1}{2}| \le \frac{1}{4}) \\ 0 & (\frac{1}{4} \le |x -\frac{1}{2}| \le \frac{1}{2}) \end{cases} $$

imax = 20
dx = 1/imax
x = np.linspace(0, 1, imax + 1)

# step 2
u = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])

# step 3
dt = 1e-3
kap = dt/dx**2
# nmax = 1000
nmax = 70
fname = "data/course1.dat"
if os.path.exists(fname):
    os.remove(fname)
f = open(fname, "a")  # 書き込み用のファイルをappendモードで開く

for n in range(nmax):
    if n % 10 == 0:
        plt.plot(x, u, label=f"t = {dt * n:.2f}")
    np.savetxt(f, u.reshape(1, len(u)))  # u[0]を末尾に書き込む
    u = fdm(u)
plt.plot(x, u, label=f"t = {dt * (n+1):.2f}")
plt.legend()
plt.xlabel("$x$")
plt.ylabel(r"$u$")
f.close()

A.2 　課程2：　 1階線形波動方程式の数値解放

例題

方程式と初期、境界条件: $$ \frac{∂u}{∂t} = - \frac{∂u}{∂x}\\ u(0,t) = u(1,t) = 0 (0<t) \\ u(x,0) = \Lambda(x) \equiv \begin{cases} 1 - 4|x -\frac{1}{2}| & (|x - \frac{1}{2}| \le \frac{1}{4}) \\ 0 & (\frac{1}{4} \le |x -\frac{1}{2}| \le \frac{1}{2}) \end{cases} $$

解析解は$u(x,t) = \Lambda(x-t) \ (t \le 0.25)$である。 ラックス・ヴェンドロフ法によって差分方程式を作ると次のようになる。

$$ \begin{aligned} u_{i}^{n+1} &= u_i^n - \frac{s}{2h}(u_{i+1}^n - u_{i-1}^n) + \frac{s^2}{2h^2}(u_{i+1}^n - 2u_i^n + u_{i-1}^n) \\ u_i^n &= u_I^n = 0 \hspace{4em} (n = 0, 1, 2, \dots) \\ u_i^0 &= \begin{cases} 1 - 4|x_i -\frac{1}{2}| & (|x_i - \frac{1}{2}| \le \frac{1}{4}) \\ 0 & (\frac{1}{4} \le |x_i -\frac{1}{2}| \le \frac{1}{2}) \end{cases} \end{aligned} $$

先の課程の最後に作成したプログラムの差分方程式を書き換える。

imax = 50
dx = 1/imax
dt = 1e-2
kap = dt/dx
nmax = 10
fname = "data/course2.dat"

x = np.linspace(0, 1, imax + 1)

# step 2
u = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
plt.plot(x + 0.1, u, label='解析解')

# テキストファイルが存在していたら削除する
if os.path.exists(fname):
    os.remove(fname)

f = open(fname, "a")  # 書き込み用のファイルをappendモードで開く

# step 3
for n in range(nmax):
    np.savetxt(f, u.reshape(1, len(u)))  # u[0]を末尾に書き込む
    if n % 10 == 0:
        plt.plot(x, u, label=f"t = {dt * n:.2f}")
    u0 = u.copy()
    for i in range(imax):
        u[i] = u0[i] - 0.5 * kap * (u0[i+1] - u0[i-1]) + 0.5 * kap**2 * (u0[i+1] - 2*u0[i] + u0[i-1])
    u[0] = 0
    u[-1] = 0

plt.plot(x, u, label=f"t = {dt * (n+1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.legend()
f.close()

ラックス・ヴェンドロフ法の関数化

コードの可読性を高め、さらに改変が容易になるようにまとまった処理をするコードを関数としてひとまとめにしておくことは、いかなるプログラムにおいても有用である。 今回のコードではラックス・ヴェンドロフ法による計算で各ステップごとに同じ処理を行っているため、これを関数としてひとまとめにしておこう。

念の為境界条件も引数として準備しておくが、境界条件を変える必要がないものは特に設定する必要はない。 もしくはグローバル変数を用いて設定しても良いだろう。

# Lax-Wendroff method
def lax_wendroff(u, bd0=0, bd1=0):
    w = np.empty(len(u))
    w[1:-1] = u[1:-1] - 0.5 * kap * (u[2:] - u[:-2]) + 0.5 * kap ** 2 * (u[2:] - 2 * u[1:-1] + u[:-2])
    w[0] = bd0
    w[-1] = bd1
    return w

imax = 50
dx = 1/imax
dt = 1e-2
kap = dt/dx
nmax = 10
fname = "data/course2v2.dat"

x = np.linspace(0, 1, imax + 1)

# step 2
u = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
plt.plot(x + 0.1, u, label='解析解')

# テキストファイルが存在していたら削除する
if os.path.exists(fname):
    os.remove(fname)

f = open(fname, "a")  # 書き込み用のファイルをappendモードで開く

# step 3
for n in range(nmax):
    np.savetxt(f, u.reshape(1, len(u)))  # u[0]を末尾に書き込む
    if (n) % 10 == 0:
        plt.plot(x, u, label=f"t = {dt * n:.2f}")
    u = lax_wendroff(u)
plt.plot(x, u, label=f"t = {dt * (n+1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.legend()
f.close()

先のコードで確認を忘れていたが、for文も関数も結局同じであることを確認しておく。

!diff -s ./data/course2v2.dat ./data/course2.dat

Files ./data/course2v2.dat and ./data/course2.dat are identical

数値不安定性

dt = 0.03

imax = 50
dt = 3e-2
kap = dt/dx
nmax = 10
dx = 1/imax

x = np.linspace(0, 1, imax + 1)

# step 2
u = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
plt.plot(x + 0.1, u, label='解析解')


# step 3
for n in range(nmax):
    if (n) % 10 == 0:
        plt.plot(x, u, label=f"t = {dt * n:.2f}")
    u = lax_wendroff(u)
plt.plot(x, u, label=f"t = {dt * (n+1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.legend()
# plt.ylim(-1, 1)
plt.plot([0.5,0.5])

[<matplotlib.lines.Line2D at 0x12e6d41c0>]

dt = 0.02

右にずれていく

imax = 50
dt = 2e-2
kap = dt/dx
nmax = 10
dx = 1/imax

x = np.linspace(0, 1, imax + 1)

# step 2
u = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])


# step 3
for n in range(nmax):
    if (n) % 5 == 0:
        plt.plot(x, u, label=f"t = {dt * n:.2f}")
    u = lax_wendroff(u)
plt.plot(x, u, label=f"t = {dt * (n+1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.legend()
plt.plot([0.5,0.5])

[<matplotlib.lines.Line2D at 0x11f97bca0>]

矩形型初期条件

$$ u(x, 0) = \Omega(x) = \begin{cases} 1 & \left( \frac{1}{3} \le x \le \frac{2}{3} \right) \\ 0 & (\text{others}) \end{cases} $$

imax = 50
dt = 1e-2
dx = 1/imax
kap = dt/dx
nmax = 10

x = np.linspace(0, 1, imax + 1)

# step 2
u = np.array([1 if np.abs(x[i] - 0.5) <= 1./6. else 0 for i in range(imax+1)])
plt.plot(x, u, label=f"t = 0")
plt.plot(x + 0.1, u, label='解析解')

# step 3
for n in range(nmax):
    u = lax_wendroff(u)
plt.plot(x, u, label=f"t = {dt * (n+1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.legend()

<matplotlib.legend.Legend at 0x12f3c2ac0>

波形の両端で誤差が生じていることがわかる。

A.3 　課程3：　　2階線型波動方程式の数値解法

例題

いよいよ2階線形波動方程式$\frac{∂^2u}{∂t^2} = \frac{∂^2 u}{∂x^2}$を数値的に解いてゆく。数学モデルは次のように、2つの一階偏微分方程式に分けられる。

$$ \begin{aligned} \frac{∂u}{∂t} &= -\frac{∂v}{∂x}, \hspace{4em} \frac{∂u}{∂t} = -\frac{∂u}{∂x} \\ u(0,t) &= u(1, t) = 0, \hspace{4em} v(0,t) = v(1, t) = 0 \hspace{4em} (0 < t < \infty) \\ u(x,0) &= v(x, 0) = \Lambda(x) \equiv \begin{cases} 1 - 4|x -\frac{1}{2}| & (|x - \frac{1}{2}| \le \frac{1}{4}) \\ 0 & (\frac{1}{4} \le |x -\frac{1}{2}| \le \frac{1}{2}) \end{cases} \end{aligned} $$

ここでベクトル $$ \boldsymbol{u} = \begin{pmatrix} u \\ v \end{pmatrix} \hspace{4em} \boldsymbol{w} = \begin{pmatrix} u \\ v \end{pmatrix} $$ を用いると次のように書き改められる：

$$ \begin{aligned} \frac{∂\boldsymbol{u}}{∂t} &= - \frac{∂\boldsymbol{w}}{∂x} \hspace{4em} (0 < x < 1, 0 < t < \infty) \\ \begin{pmatrix} u(0,t) \\ \frac{∂v}{∂x}(0, t) \end{pmatrix} &= \begin{pmatrix} u(1,t) \\ \frac{∂v}{∂x}(1, t) \end{pmatrix} = \boldsymbol{0} \hspace{4em } (0 < t < \infty)\\ \boldsymbol{u}(x, 0) &= \begin{pmatrix} \Lambda(x) \\ \Lambda(x) \end{pmatrix} \hspace{4em} (0 \le x \le 1) \end{aligned} $$

差分方程式は次のようになる:

$$ \begin{aligned} \boldsymbol{u}_{i}^{n+1} &= \boldsymbol{u}_i^n - \frac{\kappa}{2}(\boldsymbol{w}_{i+1}^n - \boldsymbol{w}_{i-1}^n) + \frac{\kappa^2}{2}(\boldsymbol{u}_{i+1}^{n} - 2\boldsymbol{u}_{i}^{n} + \boldsymbol{u}_{i-1}^{n}) \hspace{4em} (n = 0, 1, 2, \dots, i = 1, 2, \dots, I-1) \\ u_{0}^{n} &= u_{I}^{n} = 0 \\ v_0^n &= v_1^n, \hspace{4em} v_I^n = v_{I-1}^n \hspace{4em} (n = 1, 2, \dots) \\ \boldsymbol{u}_i^0 &= \begin{pmatrix} \Lambda(x_i) \\ \Lambda(x_i) \end{pmatrix} \hspace{4em} (i = 0, 1, 2, \dots, I) \end{aligned} $$

ここで差分方程式を成分別に書くと

$$ \begin{aligned} u_{i}^{n+1} &= u_i^n - \frac{\kappa}{2}(v_{i+1}^n - v_{i-1}^n) + \frac{\kappa^2}{2}(u_{i+1}^{n} - 2u_{i}^{n} + u_{i-1}^{n}) \\ v_{i}^{n+1} &= v_i^n - \frac{\kappa}{2}(u_{i+1}^n - u_{i-1}^n) + \frac{\kappa^2}{2}(v_{i+1}^{n} - 2v_{i}^{n} + v_{i-1}^{n}) \end{aligned} $$

まずは簡単な書き換え例

先のコードを二次元用に書き換える

imax = 50
nmax = 10
qmax = 2
dx = 1/imax
dt = 1e-2
kap = dt/dx
fname = "data/course3-1.dat"
n = 0

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0]
plt.plot(x, u[0], label=f"t = {dt * n:.2f}")
plt.plot(x + 0.1, u[0], label='解析解')

# テキストファイルが存在していたら削除する
if os.path.exists(fname):
    os.remove(fname)

f = open(fname, "a")  # 書き込み用のファイルをappendモードで開く
for n in range(nmax):
    np.savetxt(f, u[0].reshape(1, len(u[0])))  # u[0]を末尾に書き込む
    u0 = u.copy()
    for i in range(1, imax):
        for q in range(qmax):
            u[q][i] = np.array([u0[q-1][i] - 0.5 * kap * (u0[q-1][i+1] - u0[q-1][i-1]) + 0.5 * kap**2 * (u0[q][i+1] - 2*u0[q][i] + u0[q][i-1])])
    for q in range(qmax):
        u[q][0] = 0
        u[q][-1] = 0
    if not np.all(u[0] == u[1]):
        print('一つ目と二つ目で値が異なります')
plt.plot(x, u[0], label=f"t = {dt * n:.2f}")
plt.legend()
f.close()

文字の出力がないので、各ステップにおいてu[0]とu[1]が等しいことがわかる。

!diff -s ./data/course2.dat ./data/course3-1.dat

Files ./data/course2.dat and ./data/course3-1.dat are identical

結果は課程2、3で一致している

二次元配列に対するLax-Wendroff法の関数化

# Lax-Wendroff method　for 2-d array
def lax_w_2d(u):
    v = np.empty_like(u)
    for q in range(qmax):
        v[q][1:-1] = u[q][1:-1] - 0.5 * kap * (u[q-1][2:] - u[q-1][:-2]) + 0.5 * kap ** 2 * (u[q][2:] - 2 * u[q][1:-1] + u[q][:-2])
        v[q][0] = v[q][-1] = 0
    return v

imax = 50
nmax = 10
qmax = 2
dx = 1/imax
dt = 1e-2
kap = dt/dx
fname = "data/course3-1_v2.dat"

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()
plt.plot(x, u[0], label=f"t = 0")
plt.plot(x + 0.1, u[0], label='解析解')

# テキストファイルが存在していたら削除する
if os.path.exists(fname):
    os.remove(fname)

f = open(fname, "a")  # 書き込み用のファイルをappendモードで開く
for n in range(nmax):
    np.savetxt(f, u[0].reshape(1, len(u[0])))  # u[0]を末尾に書き込む
    u = lax_w_2d(u)
    if not np.all(u[0] == u[1]):
        print('一つ目と二つ目で値が異なります')
plt.plot(x, u[0], label=f"t = {dt * (n+1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u[0]$")
plt.legend()
f.close()

!diff -s ./data/course3-1_v2.dat ./data/course3-1.dat

Files ./data/course3-1_v2.dat and ./data/course3-1.dat are identical

続いて波動方程式を解く

(これって上で書いたコードと同じになるけど、手順分かていないのでは？)

imax = 50
nmax = 10
qmax = 2
dt = 1e-2
kap = dt/dx
dx = 1/imax
fname = "data/course3-2.dat"

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()
plt.plot(x, u[0], label=f"t = 0")
plt.plot(x + 0.1, u[0], label='解析解')

# テキストファイルが存在していたら削除する
if os.path.exists(fname):
    os.remove(fname)

f = open(fname, "a")  # 書き込み用のファイルをappendモードで開く
for n in range(nmax):
    np.savetxt(f, u[0].reshape(1, len(u[0])))  # u[0]を末尾に書き込む
    u = lax_w_2d(u)
plt.plot(x, u[0], label=f"u[0](t = {dt * (n+1):.2f})")
plt.plot(x, u[1], label=f"u[1](t = {dt * (n+1):.2f})")
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.legend()
f.close()

試してみよう1 初期条件と波の伝播方向

一つ目

imax = 50
nmax = 10
qmax = 2
dt = 1e-2
kap = dt/dx
dx = 1/imax
fname = "data/course3-3.dat"

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = -u[0].copy()
plt.plot(x, u[0], label=f"t = 0")

# テキストファイルが存在していたら削除する
if os.path.exists(fname):
    os.remove(fname)

f = open(fname, "a")  # 書き込み用のファイルをappendモードで開く
for n in range(nmax):
    np.savetxt(f, u[0].reshape(1, len(u[0])))  # u[0]を末尾に書き込む
    u = lax_w_2d(u)
plt.plot(x, u[0], label=f"u[0](t = {dt * (n+1):.2f})")
plt.plot(x, u[1], label=f"u[1](t = {dt * (n+1):.2f})")
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.legend()
f.close()

二つ目

imax = 50
nmax = 30
qmax = 2
dt = 1e-2
kap = dt/dx
dx = 1/imax
fname = "data/course3-4.dat"

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
plt.plot(x, u[0], label=f"u[0](t = 0)")
plt.plot(x, u[1], label=f"u[1](t = 0)")
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
    os.remove(fname)

# f = open(fname, "a")  # 書き込み用のファイルをappendモードで開く
for n in range(nmax):
    u = lax_w_2d(u)
    if (n+1) % 10 == 0:
        plt.plot(x, u[0], label=f"u[0](t = {dt * (n+1):.1f})")
        plt.plot(x, u[1], label=f"u[1](t = {dt * (n+1):.1f})")
    # np.savetxt(f, u[0].reshape(1, len(u[0])))  # u[0]を末尾に書き込む
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.legend()
# f.close()

<matplotlib.legend.Legend at 0x128a7cf70>

試してみよう2　　数値不安定性

imax = 50
nmax = 10
qmax = 2
dt = 3e-2
kap = dt/dx
dx = 1/imax

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()
plt.plot(x, u[0], label=f"t = 0")

for n in range(nmax):
    u = lax_w_2d(u)
plt.plot(x, u[0], label=f"u[0](t = {dt * (n + 1):.2f})")
plt.xlabel("$x$")
plt.ylabel(r"$u[0]$")
plt.legend()

<matplotlib.legend.Legend at 0x128c45cd0>

試してみよう3　 反射

# Lax-Wendroff method　for 2-d array 
def lax_w_ref(u):
    v = np.empty_like(u)
    for q in range(qmax):
        v[q][1:-1] = u[q][1:-1] - 0.5 * kap * (u[q-1][2:] - u[q-1][:-2]) + 0.5 * kap ** 2 * (u[q][2:] - 2 * u[q][1:-1] + u[q][:-2])
    v[0][0] = v[0][-1] = 0
    v[1][0] = v[1][1]
    v[1][-1] = v[1][-2]
    return v

imax = 50
nmax = 300
qmax = 2
dx = 1/imax
dt = 1e-2
kap = dt/dx

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()


fig, ax = plt.subplots(figsize=(6,4))
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0], label="t = 0")
ims = []
ax.set_xlabel("$x$")
ax.set_ylabel(r"$u[0]$")
ax2.set_xlabel("$x$")
ax2.set_ylabel(r"$u[0]$")

for n in range(nmax):
    u = lax_w_ref(u)
    txt = ax.text(0, 0.95, f"t = {dt * (n + 1):.1f}")
    im, = ax.plot(x, u[0], "black")
    ims.append([im, txt])
    if (n+1) % 30 == 0:
        ax2.plot(x, u[0], label=f"t = {dt * (n+1):.2f}")
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
ax2.legend()
plt.close(fig)                # avoid plotting a spare static plot

A.4 　課程4：　　2段階ラックス・ヴェンドロフ法

これまでは線形波動方程式をラックス・ヴェンドロフ法を用いて計算してきた。これからは非線形の計算を行うために、2段階ラックス・ヴェンドロフ法を用いたコードに書き換えていく。 2段階ラックス・ヴェンドロフ法にするには差分方程式を $$ \begin{aligned} \boldsymbol{u}_i^{\overline{n+1}} &= \boldsymbol{u}_i^n - \frac{s}{h}(\boldsymbol{w}_{i+1}^n - \boldsymbol{w}_i^n) \\ \boldsymbol{u}_i^{n+1} &= \frac{1}{2}\left[\boldsymbol{u}_i^n + \boldsymbol{u}_i^{\overline{n+1}} - \frac{s}{h}(\boldsymbol{w}_{i}^{\overline{n+1}} - \boldsymbol{w}_{i-1}^{\overline{n+1}})\right] \end{aligned} $$ とすればよい。 ここで新たな変数$\boldsymbol{w}_1^n$と$\boldsymbol{u}_i^{\overline{n+1}}$が必要となる。

手順1 　配列wの導入

imax = 50
nmax = 10
qmax = 2
dx = 1/imax
dt = 1e-2
kap = dt/dx
fname = "data/course4-1.dat"

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))

u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()

plt.plot(x, u[0], label=f"t = 0")
plt.plot(x + 0.1, u[0], label='解析解')

# テキストファイルが存在していたら削除する
if os.path.exists(fname):
    os.remove(fname)

f = open(fname, "a")  # 書き込み用のファイルをappendモードで開く
for n in range(nmax):
    np.savetxt(f, u[0].reshape(1, len(u[0])))  # u[0]を末尾に書き込む
    u0 = u.copy()
    w[1] = u[0].copy()
    w[0] = u[1].copy()
    for q in range(qmax):
        u[q][0] = 0
        u[q][-1] = 0
        for i in range(1, imax):
            u[q][i] = np.array([u0[q][i] - 0.5 * kap * (w[q][i+1] - w[q][i-1]) + 0.5 * kap**2 * (u0[q][i+1] - 2*u0[q][i] + u0[q][i-1])])

plt.plot(x, u[0], label=f"t = {dt * (n+1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u[0]$")
plt.legend()
f.close()

ここでも結果に変化がないことを確認しておく

! diff -s ./data/course2.dat ./data/course4-1.dat

Files ./data/course2.dat and ./data/course4-1.dat are identical

手順2 　2段階ラックス・ヴェンドロフ法の差分式への変更

imax = 50
nmax = 10
qmax = 2
dx = 1/imax
dt = 1e-2
kap = dt/dx
fname = "data/course4-2.dat"

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))

u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()

plt.plot(x, u[0], label=f"t = 0")
plt.plot(x + 0.1, u[0], label='解析解')

# テキストファイルが存在していたら削除する
if os.path.exists(fname):
    os.remove(fname)

f = open(fname, "a")  # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
    np.savetxt(f, u[0].reshape(1, len(u[0])))  # u[0]を末尾に書き込む
    u0 = u.copy()
    w[1] = u[0].copy()
    w[0] = u[1].copy()
    for i in range(0, imax):
        for q in range(qmax):
            uh[q][i] = u0[q][i] - kap * (w[q][i+1] - w[q][i])
    uh[0][imax] = 0.0
    uh[1][imax] = uh[1][imax - 1]
    w[0] = uh[1].copy()
    w[1] = uh[0].copy()
    for i in range(1, imax + 1):
        for q in range(qmax):
            u[q][i] = 0.5 * (u0[q][i] + uh[q][i] - kap * (w[q][i] - w[q][i-1]))
    u[0][0] = 0.0
    u[1][0] = u[1][1]

plt.plot(x, u[1], label=f"t = {dt * (n+1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.legend()
f.close()

2段階ラックス・ヴェンドロフ法の関数化

今回も積分の計算は関数としてひとまとめにしておこう。 関数化するにあたって、wの計算が必要になるため、それをcalw()として引数にとることにする。

また、先の関数では境界条件を引数に設定していたが、ここでは手順4で課す周期的境界条件を課しておくことにする。

# Lax-Wendroff two step method
def lax_w_2step(u, calw):
    # Firts step
    w = calw(u)
    u_first = np.empty_like(u)
    u_first[:,0:-1] = u[:,0:-1] - kap * (w[:, 1:] - w[:,0:-1]) 
    u_first[:, -1] = u_first[:, 0] # BdC
    # Second step
    w = calw(u_first)
    u_second = np.empty_like(u)
    u_second[:, 1:] = 0.5 * (u[:, 1:] + u_first[:, 1:] - kap * (w[:, 1:] - w[:, :-1]))
    u_second[:, 0] = u_second[:, -1] # BdC
    return u_second

手順3　　$w$の計算の関数化

実習資料に合わせるならば

def calw(u, w):
    w[0] = u[1].copy()
    w[1] = u[0].copy()

とすれば良いだろう。

しかし、このような関数内部で変数を変更しようとすると、下の「おまけ」で見るように代入するもの（配列、文字列、配列の要素、etc.）によって処理が変わってしまうので、仕様を理解していないと意図しない挙動になりそうである。 （配列の要素の代入はきちんと置き換えられるため上の関数でも良いのだが）

そこで、次のように0、1要素を入れ替えた配列(ndarray)を返す関数を作り、wに代入するほうが良いのではないだろうか。

def calw_wave(u):
    w = np.empty_like(u)
    w[0] = u[1].copy()
    w[1] = u[0].copy()
    return w
# example of usage
# w = calw_wave(u)

この関数の実行例を見て見ると、きちんと第０要素と第１要素が入れ替わっていることがわかる。

# Test
u = np.arange(0,10).reshape(2, 5)
w = calw_wave(u)
print(f"u = \n{u}")
print(f"w = \n{w}")
print(f"u[0] == w[1]: {np.all(u[0] == w[1])}")
print(f"u[1] == w[0]: {np.all(u[1] == w[0])}")

u = 
[[0 1 2 3 4]
 [5 6 7 8 9]]
w = 
[[5 6 7 8 9]
 [0 1 2 3 4]]
u[0] == w[1]: True
u[1] == w[0]: True

おまけ　　関数内での引数の変更のテスト

ndarrayの要素に代入（変更される）

def f(x,y):
    for i in range(10):
        y[i] = x[i]
    print(f"関数内：{y}")
x = np.arange(0,10,1)
y = np.zeros(10)
print(f"初め：{y}")
f(x, y)
print(f"関数実行後：{y}")

初め：[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
関数内：[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
関数実行後：[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]

ndarray全体を代入（変更されない）

def f(x,y):
    y = x.copy()
    print(f"関数内：{y}")
x = np.arange(0,10,1)
y = np.zeros(10)
print(f"初め：{y}")
f(x, y)
print(f"関数実行後：{y}")

初め：[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
関数内：[0 1 2 3 4 5 6 7 8 9]
関数実行後：[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]

文字列（変更されない）

def f(a):
    a += "def"
    print(f"関数内：{a}")
a = "abc"
print(f"初め：{a}")
f(a)
print(f"関数実行後：{a}")

初め：abc
関数内：abcdef
関数実行後：abc

多次元のndarrayの要素の配列を代入（変更される）

def f(x, y):
    y[0] = x[1].copy()
    y[1] = x[0].copy()
x = np.arange(0,10).reshape(2,5)
y = np.zeros(10).reshape(2,5)
print(f"初め： y =\n{y}")
f(x, y)
print(f"関数実行後： y = \n{y}")

初め： y =
[[0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]]
関数実行後： y = 
[[5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4.]]

定義した関数の利用

imax = 50
nmax = 10
qmax = 2
dx = 1/imax
dt = 1e-2
kap = dt/dx
fname = "data/course4-2_v2.dat"

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))

u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()

plt.plot(x, u[0], label=f" t = 0")
plt.plot(x + 0.1, u[0], label='解析解')

# テキストファイルが存在していたら削除する
if os.path.exists(fname):
    os.remove(fname)

f = open(fname, "a")  # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
    np.savetxt(f, u[0].reshape(1, len(u[0])))  # u[0]を末尾に書き込む
    u = lax_w_2step(u, calw_wave)

plt.plot(x, u[0], label=f"t = {dt * (n+1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u[0]$")
plt.legend()
f.close()

!diff -s ./data/course4-2_v2.dat ./data/course4-2.dat

Files ./data/course4-2_v2.dat and ./data/course4-2.dat are identical

今の計算範囲では境界条件が使われていないため、関数版(*_v2.dat)とfor文版で結果が同じであることがわかる。

手順4　　周期境界条件

ステップ数を増やす nmax = 50

imax = 50
nmax = 50
qmax = 2
dx = 1/imax
dt = 1e-2
kap = dt/dx

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))

u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()

plt.plot(x, u[0], label=f"t = 0")

for n in range(nmax):
    u0 = u.copy()
    w[1] = u[0].copy()
    w[0] = u[1].copy()
    for i in range(0, imax):
        for q in range(qmax):
            uh[q][i] = u0[q][i] - kap * (w[q][i+1] - w[q][i])
    uh[0][imax] = 0.0
    uh[1][imax] = uh[1][imax - 1]
    w[0] = uh[1].copy()
    w[1] = uh[0].copy()
    for i in range(1, imax + 1):
        for q in range(qmax):
            u[q][i] = 0.5 * (u0[q][i] + uh[q][i]) - kap * (w[q][i] - w[q][i-1])
    u[0][0] = 0.0
    u[1][0] = u[1][1]
    
plt.plot(x, u[0], label=f"t = {dt * (n + 1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u[0]$")
plt.legend()

<matplotlib.legend.Legend at 0x12f971880>

関数の書き換えの手間を省くため、ここでは関数化前のコードを利用した。
形が崩れていることがわかる。

周期的境界条件の導入

一周するようにnmaxを100に変えておく。

imax = 50
nmax = 100
qmax = 2
dx = 1/imax
dt = 1e-2
kap = dt/dx
fname = "data/course4-4.dat"

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))

u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()
# plt.plot(x, u[0], label=f"t = 0")


fig, ax = plt.subplots(figsize=(6,4))
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0], label="t = 0")
ims = []
ax.set_xlabel("$x$")
ax.set_ylabel(r"$u[0]$")
ax2.set_xlabel("$x$")
ax2.set_ylabel(r"$u[0]$")


# テキストファイルが存在していたら削除する
if os.path.exists(fname):
    os.remove(fname)
f = open(fname, "a")  # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
    np.savetxt(f, u[0].reshape(1, len(u[0])))  # u[0]を末尾に書き込む
    u = lax_w_2step(u, calw_wave)
    txt = ax.text(0, 0.95, f"t = {dt * (n + 1):.1f}")
    im, = ax.plot(x, u[0], "black")
    ims.append([im, txt])
    if (n+1) % 30 == 0:
        ax2.plot(x, u[0], label=f"t = {dt * (n+1):.2f}")
f.close()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
ax2.legend()
plt.close(fig)                # avoid plotting a spare static plot

アニメーションから、右端に達した波が左端から出てきていることがよくわかる。

A.5 課程5:　　非線形方程式<流体力学方程式>の数値解法

いよいよ非線形方程式の数値解法に入る。

手順1　 変数$\boldsymbol{u}$の数qmaxを3に増やす

まず、qmaxの値を3に増やし、増やした2番目の要素は1番目の要素と同じになるように書き換える。 つまり

u[2][i] = u[1][i]

となるようにする。

この時、先に定義したcalw_wave()では対応できないため、これも少し修正する。(もっとも、ここでしか使わないため適当である。)

def calw_wave_q3(u):
    w = np.empty_like(u)
    w[0] = u[1].copy()
    w[1] = u[0].copy()
    w[2] = w[1].copy()
    return w

imax = 50
nmax = 100
qmax = 3
dx = 1/imax
dt = 1e-2
kap = dt/dx
fname = "data/course5-1.dat"

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))

u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()

# テキストファイルが存在していたら削除する
if os.path.exists(fname):
    os.remove(fname)

f = open(fname, "a")  # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
    if (n) % 20 == 0:
        plt.plot(x, u[0], label=f"t = {dt * n:.2f}")
    np.savetxt(f, u[0].reshape(1, len(u[0])))  # u[0]を末尾に書き込む
    u = lax_w_2step(u, calw_wave_q3)
        
plt.plot(x, u[0], label=f"t = {dt * (n+1):.1f}")
plt.xlabel("$x$")
plt.ylabel(r"$u[0]$")
plt.legend(loc="upper left")
f.close()

!diff -s ./data/course4-4.dat ./data/course5-1.dat

Files ./data/course4-4.dat and ./data/course5-1.dat are identical

計算結果が先のものと変わらないことも確かめられた。

手順2 $\boldsymbol{w}$を流体力学方程式の流速密度にする

流速密度行列$\boldsymbol{w}$の計算を $$ \boldsymbol{w} = \begin{pmatrix} \rho v \\ \rho v^2 + p \\ (e + p) v \end{pmatrix} $$ に対応するように書き換える。ここで、$\rho$は気体の密度、$v$は気体の速度、$e$はエネルギー密度である。 この時、各メッシュにおいて $$ \boldsymbol{u} = \begin{pmatrix} \rho \\ \rho v \\ e \end{pmatrix} $$ の成分を用いて $$ u = \frac{\rho v}{\rho}, \hspace{4em} p = (\gamma - 1)\left(e - \frac{\rho}{2}v^2 \right) $$ と計算する。

def calw_fluid(u):
    w = np.empty_like(u)
    rho = u[0].copy()
    rhov = u[1].copy()
    e = u[2].copy()
    v = rhov / rho
    p = (gamma - 1) * (e - rhov * v / 2)
    w[0] = rhov
    w[1]= rhov * v + p
    w[2] = (e + p) * v
    return w

手順3 初期条件の設定

$$ \begin{aligned} \begin{array}{ll} \rho=1 & (0 \leq x \leq 1) \\ \rho v=0 & (0 \leq x \leq 1) \\ e=\left\{\begin{array}{cc} 1.9-4\left|x-\frac{1}{2}\right| & \left(\frac{1}{4} \leq x \leq \frac{3}{4}\right) \\ 0.9 & \left(0 \leq x<\frac{1}{4}, \frac{3}{4}<x \leq 1\right) \end{array}\right. \end{array} \end{aligned} $$

imax = 50
nmax = 50
qmax = 3
dx = 1/imax
dt = 1e-2
kap = dt/dx
gamma = 5/3

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))

u[0] = np.ones(imax+1)
u[1] = np.zeros(imax+1)
u[2] = np.array([1.9 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0.9 for i in range(imax + 1)])

fig, axes = plt.subplots(1,3,figsize=(12, 3))  # for statical plot
fig_anim, axes_anim = plt.subplots(1,3,figsize=(15, 4))  # for animation plot
ims = []
plots = [0] * 3

for i, ax in enumerate(axes):
    ax.plot(x, u[i], label="t = 0")

for n in range(nmax):
    # animation
    for i, ax_anim in enumerate(axes_anim):
        plots[i] = ax_anim.plot(x, u[i], c="black")
    txt = axes_anim[0].text(0, 1.5, f"t = {dt * (n):.1f}", size=16)
    plots[0] = [plots[0][0],txt]  # to add text in the first ax, this is needed
    ims.append(plots[0]+plots[1]+plots[2])  # for subplots in animation
    
    u = lax_w_2step(u, calw_fluid)
    # nomal plot
    if (n + 1) % 10 == 0:
        for i, ax in enumerate(axes):
            ax.plot(x, u[i], label=f"t = {dt * (n+1):.1f}")
# Set x, and y labels
for axs in [axes, axes_anim]:
    axs[1].set_xlabel("$x$")
    axs[0].set_ylabel(r"$\rho$")
    axs[1].set_ylabel(r"$\rho v$")
    axs[2].set_ylabel(r"energy")
axes[0].legend(loc="best")
# Adjust layouts
fig.tight_layout()
fig_anim.tight_layout()
# Start animation
anim = ArtistAnimation(fig_anim, ims, interval=100)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
plt.close(fig_anim)    # close animation plot  

試してみよう1 音波のシミュレーション

A = 0.001

imax = 50
nmax = 100
qmax = 3
dx = 1/imax
dt = 1e-2
kap = dt/dx
gamma = 5/3
A = 1e-3

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))

u[0] = 1 + A * np.cos(2*np.pi*x)
u[1] = A * np.cos(2*np.pi*x)
u[2] = 0.9 * (1 + gamma * A * np.cos(2*np.pi*x))
# plt.plot(x, u[0], label=f"t = 0")


fig, ax = plt.subplots(figsize=(6,4))
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0], label="t = 0")
ims = []
ax.set_xlabel("$x$")
ax.set_ylabel(r"$u[0]$")
ax2.set_xlabel("$x$")
ax2.set_ylabel(r"$u[0]$")


# テキストファイルが存在していたら削除する
if os.path.exists(fname):
    os.remove(fname)
f = open(fname, "a")  # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
    np.savetxt(f, u[0].reshape(1, len(u[0])))  # u[0]を末尾に書き込む
    u = lax_w_2step(u, calw_fluid) # 積分1ステップ進める
    # animation
    txt = ax.text(0, 0.95, f"t = {dt * (n + 1):.1f}")
    im, = ax.plot(x, u[0], "black")
    ims.append([im, txt])
    if (n+1) % 25 == 0:
        ax2.plot(x, u[0], ".", label=f"t = {dt * (n+1):.2f}")
f.close()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
ax2.legend()
plt.close(fig)   

この時には安定した波動が観察される。

A = 0.1

この時は非線形の影響が強く出ることになる。

imax = 50
nmax = 200
qmax = 3
dx = 1/imax
dt = 1e-2
kap = dt/dx
gamma = 5/3
A = 0.1
fname = "data/course5-ex1-A01.dat"

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))

u[0] = 1 + A * np.sin(2*np.pi*x)
u[1] = A * np.sin(2*np.pi*x)
u[2] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*x))
# plt.plot(x, u[0], label=f"t = 0")


fig, ax = plt.subplots(figsize=(6,4))
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0], label="t = 0")
ims = []
ax.set_xlabel("$x$")
ax.set_ylabel(r"$u[１]$")
ax2.set_xlabel("$x$")
ax2.set_ylabel(r"$u[1]$")


# テキストファイルが存在していたら削除する
if os.path.exists(fname):
    os.remove(fname)
f = open(fname, "a")  # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
    np.savetxt(f, u[0].reshape(1, len(u[1])))  # u[0]を末尾に書き込む
    u = lax_w_2step(u, calw_fluid) # 積分1ステップ進める
    # animation
    txt = ax.text(0, 0.3, f"t = {dt * (n + 1):.1f}")
    im, = ax.plot(x, u[0], "black")
    ims.append([im, txt])
    if (n+1) % 100 == 0:
        ax2.plot(x, u[0], label=f"t = {dt * (n+1):.2f}")
f.close()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
ax2.legend()
plt.close(fig)   

Font 'rm' does not have a glyph for '\uff11' [U+ff11], substituting with a dummy symbol.

方程式の非線形な部分が効いており、衝撃波のようなものが観察される。

[解決]関数が間違えている可能性を考慮してfor文を回す

初期設定が、実習書では$\cos$、教科書では$\sin$が使われていたが、それに気づかず実行結果が異なって見えたため、関数化の際の間違いがないかを確認するべくfor文を回して確認した。

実際は初期設定を変えることで解決したが、関数化の妥当性を見るために以下の記述を残しておく。

# for文を使った計算（検算用）
imax = 50
nmax = 200
qmax = 3
dx = 1/imax
dt = 1e-2
kap = dt/dx
gamma = 5/3
A = 0.1
fname = "data/course5-ex1-A01_check.dat"

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))

u[0] = 1 + A * np.sin(2*np.pi*x)
u[1] = A * np.sin(2*np.pi*x)
u[2] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*x))
# plt.plot(x, u[0], label=f"t = 0")


fig, ax = plt.subplots(figsize=(6,4))
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[1], label="t = 0")
ims = []
ax.set_xlabel("$x$")
ax.set_ylabel(r"$u[１]$")
ax2.set_xlabel("$x$")
ax2.set_ylabel(r"$u[1]$")


# テキストファイルが存在していたら削除する
if os.path.exists(fname):
    os.remove(fname)
f = open(fname, "a")  # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
    np.savetxt(f, u[1].reshape(1, len(u[1])))  # u[1]を末尾に書き込む
    u0 = u.copy()
    w = calw_fluid(u)
    for i in range(0, imax):
        for q in range(qmax):
            uh[q][i] = u0[q][i] - kap * (w[q][i+1] - w[q][i])
    uh[:, -1] = uh[:,0]
    w = calw_fluid(uh)
    for i in range(1, imax + 1):
        for q in range(qmax):
            u[q][i] = 0.5 * (u0[q][i] + uh[q][i] - kap * (w[q][i] - w[q][i-1]))
    u[:, 0] = u[:, -1]
    # animation
    txt = ax.text(0, 0.95, f"t = {dt * (n + 1):.1f}")
    im, = ax.plot(x, u[1], "black")
    ims.append([im, txt])
    if (n+1) % 100 == 0:
        ax2.plot(x, u[1], label=f"t = {dt * (n+1):.2f}")
f.close()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
ax2.legend()
plt.close(fig)  

Font 'rm' does not have a glyph for '\uff11' [U+ff11], substituting with a dummy symbol.

!diff -s ./data/course5-ex1-A01_check.dat ./data/course5-ex1-A01.dat

Files ./data/course5-ex1-A01_check.dat and ./data/course5-ex1-A01.dat are identical

試してみよう2 爆縮

初期条件

$$ \begin{aligned} \begin{array}{ll} \rho=1 & (0 \leq x \leq 1) \\ \rho v=0 & (0 \leq x \leq 1) \\ e=\left\{\begin{array}{cc} 0.4 + 2 \left|x-\frac{1}{2}\right| & \left(\frac{1}{4} \leq x \leq \frac{3}{4}\right) \\ 0.9 & \left(0 \leq x<\frac{1}{4}, \frac{3}{4}<x \leq 1\right) \end{array}\right. \end{array} \end{aligned} $$

コード

imax = 50
nmax = 100
qmax = 3
dx = 1/imax
dt = 1e-2
kap = dt/dx
gamma = 5/3

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))

u[0] = np.ones(imax+1)
u[1] = np.zeros(imax+1)
u[2] = np.array([0.4 + 2 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0.9 for i in range(imax + 1)])

fig, axes = plt.subplots(1,3,figsize=(12, 3))  # for static plot
fig_anim, axes_anim = plt.subplots(1,3,figsize=(15, 4))  # for animation plot
ims = []
plots = [0] * 3
labels = [r"$\rho$", r"$\rho v$", r"energy"]

# set first plot and labesl
for i, (ax, ax_anim) in enumerate(zip(axes, axes_anim)):
    ax.plot(x, u[i], label="t = 0")
    plots[i] = ax_anim.plot(x, u[i], c="black")
    
    ax.set_ylabel(labels[i])
    ax_anim.set_ylabel(labels[i])
    
    ax.set_xlabel(r"$x$")
    ax_anim.set_xlabel(r"$x$")
    
txt = axes_anim[0].text(0, 1.8, f"t = 0", size=16)
plots[0] = [plots[0][0],txt]  # to add text in the first ax, this is needed
ims.append(plots[0]+plots[1]+plots[2])  # for subplots in animation

for n in range(nmax):
    u = lax_w_2step(u, calw_fluid)
    
    # animation
    for i, ax_anim in enumerate(axes_anim):
        plots[i] = ax_anim.plot(x, u[i], c="black")
    txt = axes_anim[0].text(0, 1.5, f"t = {dt * (n):.1f}", size=16)
    plots[0] = [plots[0][0],txt]  # to add text in the first ax, this is needed
    ims.append(plots[0]+plots[1]+plots[2])  # for subplots in animation
    
    # nomal plot
    if (n + 1) % 50 == 0:
        for i, ax in enumerate(axes):
            ax.plot(x, u[i], label=f"t = {dt * (n+1):.1f}")

axes[1].legend(loc="best")
# Adjust layouts
fig.tight_layout()
fig_anim.tight_layout()
# Start animation
anim = ArtistAnimation(fig_anim, ims, interval=100)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
plt.close(fig_anim)    # close animation plot  

試してみよう3 高密度ガスの衝突

初期条件

$$ \begin{aligned} \rho &= \begin{cases} 1 + |\sin(6\pi x)| & \left(\frac{1}{6} < x < \frac{1}{3}, \frac{2}{3} < x < \frac{5}{6} \right) \\ 1 & (\textrm{others}) \end{cases} \\ v &= \begin{cases} -\sin(6\pi x) & \left(\frac{1}{6} < x < \frac{1}{3}, \frac{2}{3} < x < \frac{5}{6} \right) \\ 1 & (\textrm{others}) \end{cases}\\ p &= 0.9 \end{aligned} $$

コード

imax = 50
nmax = 200
qmax = 3
dx = 1/imax
dt = 1e-3
kap = dt/dx
gamma = 5/3

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
rho = np.ones(imax+1)
v = np.zeros(imax+1)
p = 0.9 * np.ones_like(x)
              
rho[round(1/6/dx):round(1/3/dx)] = 1 + np.abs(np.sin(6*np.pi*x[round(1/6/dx):round(1/3/dx)]))
rho[round(2/3/dx):round(5/6/dx)] = 1 + np.abs(np.sin(6*np.pi*x[round(2/3/dx):round(5/6/dx)]))
v[round(1/6/dx):round(1/3/dx)] = -np.sin(6*np.pi*x[round(1/6/dx):round(1/3/dx)])
v[round(2/3/dx):round(5/6/dx)] = -np.sin(6*np.pi*x[round(2/3/dx):round(5/6/dx)])
u[0] = rho
u[1] = rho * v
u[2] = p/(gamma-1) + rho * v * v / 2


fig, axes = plt.subplots(1,3,figsize=(12, 3))  # for static plot
fig_anim, axes_anim = plt.subplots(1,3,figsize=(15, 4))  # for animation plot
ims = []
plots = [0] * 3
labels = [r"$\rho$", r"$\rho v$", r"energy"]

# set first plot and labesl
for i, (ax, ax_anim) in enumerate(zip(axes, axes_anim)):
    ax.plot(x, u[i], label="t = 0")
    plots[i] = ax_anim.plot(x, u[i], c="black")
    
    ax.set_ylabel(labels[i])
    ax_anim.set_ylabel(labels[i])
    
    ax.set_xlabel(r"$x$")
    ax_anim.set_xlabel(r"$x$")
    
txt = axes_anim[0].text(0, 1.8, f"t = 0", size=16)
plots[0] = [plots[0][0],txt]  # to add text in the first ax, this is needed
ims.append(plots[0]+plots[1]+plots[2])  # for subplots in animation

for n in range(nmax):
    u = lax_w_2step(u, calw_fluid)
    
    # animation
    for i, ax_anim in enumerate(axes_anim):
        plots[i] = ax_anim.plot(x, u[i], c="black")
    txt = axes_anim[0].text(0, 1.8, f"t = {dt * (n):.1f}", size=16)
    plots[0] = [plots[0][0],txt]  # to add text in the first ax, this is needed
    ims.append(plots[0]+plots[1]+plots[2])  # for subplots in animation
    
    # nomal plot
    if (n + 1) % 100 == 0:
        for i, ax in enumerate(axes):
            ax.plot(x, u[i], label=f"t = {dt * (n+1):.1f}")

axes[1].legend(loc="best")
# Adjust layouts
fig.tight_layout()
fig_anim.tight_layout()
# Start animation
anim = ArtistAnimation(fig_anim, ims, interval=10)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
plt.close(fig_anim)    # close animation plot  

検証

結果がおかしいので、wの計算が積分がおかしいかどうかfor文を用いたものと比較してみる。

したのをやってみて思ったが、積分の計算は同じ結果になるのは当たり前だと思う。 wの計算も間違えようがないとも思うし、一体いどうして計算が破綻するのだろうか。

# Lax-Wendroff two step method　for ver.
def lax_w_2step_for(u, calw):
    # Firts step
    w = calw(u)
    u_first = np.empty_like(u)
    for i in range(imax):
        u_first[:,i] = u[:,i] - kap * (w[:,i+1] - w[:, i])
    u_first[:, -1] = u_first[:, 0] # BdC
    # Second step
    w = calw(u_first)
    u_second = np.empty_like(u)
    for i in range(1, imax+1):
        u_second[:, i] = 0.5 * (u[:, i] + u_first[:, i] - kap * (w[:, i] - w[:, i-1]))
    u_second[:, 0] = u_second[:, -1] # BdC
    return u_second

imax = 50
nmax = 350
qmax = 3
dx = 1/imax
dt = 1e-4
kap = dt/dx
gamma = 5/3

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))
rho = np.ones(imax+1)
v = np.ones(imax+1)
p = np.array([0.9]*(imax+1))
              
rho[round(1/6/dx):round(1/3/dx)] = 1 + np.abs(np.sin(6*np.pi*x[round(1/6/dx):round(1/3/dx)]))
rho[round(2/3/dx):round(5/6/dx)] = 1 + np.abs(np.sin(6*np.pi*x[round(2/3/dx):round(5/6/dx)]))
v[round(1/6/dx):round(1/3/dx)] = -np.sin(6*np.pi*x[round(1/6/dx):round(1/3/dx)])
v[round(2/3/dx):round(5/6/dx)] = -np.sin(6*np.pi*x[round(2/3/dx):round(5/6/dx)])
u[0] = rho
u[1] = rho * v
u[2] = p/(gamma-1) + rho * v * v / 2

v = u.copy()

# plt.plot(x, u[2], label=f"t = 0")

for n in range(nmax):
    u = lax_w_2step(u, calw_fluid)
    v = lax_w_2step_for(v, calw_fluid)
    if np.any(v != u):
        print("u is not equal tot\ v!!!")

付録B ２次元流体力学方程式の数値解法

方程式

今までは1次元の問題のみを扱ってきたが、これからは多次元の問題に取り組んでいく。初めは2次元である。

2次元流体力学方程式は次のように表せる。

$$ \begin{aligned} \begin{gathered} \frac{\partial \mathbf{u}}{\partial t}=-\frac{\partial \mathbf{w}_{1}}{\partial x}-\frac{\partial \mathbf{w}_{2}}{\partial y} \\ \mathbf{u} \equiv\left(\begin{array}{c} \rho \\ \rho v_{x} \\ \rho v_{y} \\ e \end{array}\right), \mathbf{w}_{1} \equiv\left(\begin{array}{c} \rho v_{x} \\ \rho v_{x}^{2}+p \\ \rho v_{x} v_{y} \\ (e+p) v_{x} \end{array}\right), \mathbf{w}_{2} \equiv\left(\begin{array}{c} \rho v_{y} \\ \rho v_{x} v_{y} \\ \rho v_{y}^{2}+p \\ (e+p) v_{y} \end{array}\right) \end{gathered} \end{aligned} $$

$$ e=\frac{\rho}{2} v^{2}+\frac{p}{\gamma-1} $$

はエネルギー密度

コードの書き換え

qmaxを4へ増やす

3番目の配列は全て2番目と同じにする。
上の音波のシミュレーションのコードを用いよう。

なお、アニメーションは邪魔なので消しておく。

def calw_fluid_q4(u):
    w = np.empty_like(u)
    rho = u[0].copy()
    rhov = u[1].copy()
    e = u[2].copy()
    v = rhov / rho
    p = (gamma - 1) * (e - rhov * v / 2)
    w[0] = rhov
    w[1]= rhov * v + p
    w[2] = (e + p) * v
    w[3] = w[2].copy()
    return w

imax = 50
nmax = 200
qmax = 4
dx = 1/imax
dt = 1e-2
kap = dt/dx
gamma = 5/3
A = 0.1
fname = "data/courseB-1-A03.dat"

x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))

u[0] = 1 + A * np.sin(2*np.pi*x)
u[1] = A * np.sin(2*np.pi*x)
u[2] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*x))
u[3] = u[2].copy()

# plt.plot(x, u[0], label=f"t = 0")


fig, ax = plt.subplots(figsize=(6,4))
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[1], label="t = 0")
ims = []
ax.set_xlabel("$x$")
ax.set_ylabel(r"$u[１]$")
ax2.set_xlabel("$x$")
ax2.set_ylabel(r"$u[1]$")


# テキストファイルが存在していたら削除する
if os.path.exists(fname):
    os.remove(fname)
f = open(fname, "a")  # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
    np.savetxt(f, u[1].reshape(1, len(u[1])))  # u[0]を末尾に書き込む
    u = lax_w_2step(u, calw_fluid_q4) # 積分1ステップ進める
    # animation
    txt = ax.text(0, 0.3, f"t = {dt * (n + 1):.1f}")
    im, = ax.plot(x, u[1], "black")
    ims.append([im, txt])
    if (n+1) % 100 == 0:
        ax2.plot(x, u[1], label=f"t = {dt * (n+1):.2f}")
f.close()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
ax2.legend()
plt.close(fig)   

Font 'rm' does not have a glyph for '\uff11' [U+ff11], substituting with a dummy symbol.

!diff -s data/course5-ex1-A03.dat data/courseB-1-A03.dat

Files data/course5-ex1-A03.dat and data/courseB-1-A03.dat are identical

今回の書き換えではu[0]に変化がない。

ラックス・ウェンドロフ法の関数を書き換える。

新たな次元が加わることで、ラックス・ウェンドロフ法の関数も書き換えなければならない。 変更は以下の要領で行う。

（第１段階，予測子）

uh[q][i][j] = uo[q][i][j] - kap1*(w1[q][i+1][j]-w1[q][i][j]) - kap2*(w2[q][i][j+1]-w2[q][i][j])

（第２段階，修正子）

u[q][i][j] = 0.5*(uo[q][i][j] + uh[q][i][j] - kap1*(w1[q][i][j]-w1[q][i-1][j]) 
                  - kap2*(w2[q][i][j]-w2[q][i][j-1]))

# Lax-Wendroff two step method for 2-dimention
def lax_w_2step_2d(u, calw):
    # Firts step
    w1, w2 = calw(u)
    u_first = np.ones_like(u)
    u_first[:,:-1, :-1] = u[:,:-1, :-1] - kap1 * (w1[:, 1:, :-1] - w1[:,:-1, :-1]) - kap2 * (w2[:, :-1, 1:] - w2[:,:-1, :-1])
    u_first[:, :, -1] = u_first[:, :, 0].copy() # BdC
    u_first[:, -1, :] = u_first[:, 0, :].copy() # BdC
    # Second step
    w1, w2 = calw(u_first)
    u_second = np.ones_like(u)
    u_second[:, 1:, 1:] = 0.5 * (u[:, 1:, 1:] + u_first[:, 1:, 1:] - kap1 * (w1[:, 1:, 1:] - w1[:, :-1, 1:]) - kap2 * (w2[:, 1:, 1:] - w2[:, 1:, :-1]))
    u_second[:, :, 0] = u_second[:, :, -1].copy() # BdC
    u_second[:, 0, :] = u_second[:, -1, :].copy() # BdC
    return u_second

新しく加わる独立変数yに対応するu，wの次元を増やす。

def calw_fluid_q4v2(u):
    w = np.empty_like(u)
    rho = u[0].copy()
    rhov = u[1].copy()
    e = u[2].copy()
    v = rhov / rho
    p = (gamma - 1) * (e - rhov * v / 2)
    w[0] = rhov
    w[1]= rhov * v + p
    w[2] = (e + p) * v
    w[3] = w[2].copy()
    return w, w

imax = 50
jmax = 50
nmax = 200
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 1e-2
kap1 = dt/dx
kap2 = dt/dy
# kap = dt/dx
gamma = 5/3
A = 0.1
fname = "data/courseB-2-A01.dat"
j_flag = 0


x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1) # y軸方向の
u = np.zeros((qmax, jmax+1, imax+1))  # u[q,y,x]  次元を増やすため、jmax+1 の部分を足した。
X, Y = np.meshgrid(x, y)


u[0] = 1 + A * np.sin(2*np.pi*X)
u[1] = A * np.sin(2*np.pi*X)
u[2] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*X))
u[3] = u[2].copy()

fig, ax = plt.subplots(figsize=(6,4))
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0,0,:], label="t = 0")
ims = []
ax.set_xlabel("$x$")
ax.set_ylabel(r"$u[0]$")
ax2.set_xlabel("$x$")
ax2.set_ylabel(r"$u[0]$")


# テキストファイルが存在していたら削除する
if os.path.exists(fname):
    os.remove(fname)
f = open(fname, "a")  # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
    np.savetxt(f, u[1, 0, :].reshape(1, len(u[1, 0, :]))) 
    u = lax_w_2step_2d(u, calw_fluid_q4v2) # 積分1ステップ進める
    if (n+1) % 100 == 0:
        ax2.plot(x, u[0, 0, :], label=f"t = {dt * (n+1):.2f}")
    # animation
    txt = ax.text(0, 1, f"t = {dt * (n + 1):.1f}")
    im, = ax.plot(x, u[0,0, :], "black")
    ims.append([im, txt])
f.close()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
ax2.legend()
plt.close(fig)  

 !diff -s data/courseB-2-A01.dat data/course5-ex1-A01.dat

3次元プロットも追加してみる

from mpl_toolkits.mplot3d import Axes3D としてmatplotlib　の3Dプロットを用いる。

しかし、描画にものすごく時間がかかるので、3Dプロットは今後使わず、等高線プロットにしておこう。

imax = 50
jmax = 50
nmax = 50
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 1e-2
kap1 = dt/dx
kap2 = dt/dy
gamma = 5/3
A = 0.1

x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1) # y軸方向の
X, Y = np.meshgrid(x, y) 
u = np.zeros((qmax, imax+1, jmax+1))  # 次元を増やすため、jmax+1 の部分を足した。

u[0] = 1 + A * np.sin(2*np.pi*X)
u[1] = A * np.sin(2*np.pi*X)
u[2] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*X))
u[3] = u[2].copy()

# fig, ax = plt.subplots(figsize=(6,4))
fig2, ax2 = plt.subplots(figsize=(6,4))
fig3, ax3 = plt.subplots(figsize=(9,6),subplot_kw=dict(projection='3d')) # for 3d animation
ax2.plot(x, u[1,0,:], label="t = 0")
ims = []
# ax.set_xlabel("$x$")
# ax.set_ylabel(r"$u[１]$")
ax2.set_xlabel("$x$")
ax2.set_ylabel(r"$u[1]$")
ax3.set_xlabel("$x$")
ax3.set_ylabel("$y$")
ax3.set_zlabel("u[1]")

for n in range(nmax):
    u = lax_w_2step_2d(u, calw_fluid_q4v2) # 積分1ステップ進める
    if (n+1) % 100 == 0:
        ax2.plot(x, u[1, 0, :], label=f"t = {dt * (n+1):.2f}")
    txt = ax3.text(0, 0, 0, f"t = {dt * (n + 1):.1f}")
    im = ax3.plot_surface(X, Y, u[0, :, :], cmap=cm.viridis)
    ims.append([im, txt])
anim = ArtistAnimation(fig3, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
ax2.legend()
plt.close(fig3)  

等高線プロットを作成してみる。

等高線プロットはax.contourで作成する。

これはこれで時間がかかるけれど、3次元プロットよりも幾分かマシである。

imax = 50
jmax = 50
nmax = 200
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 1e-2
kap1 = dt/dx
kap2 = dt/dy
gamma = 5/3
A = 0.1

x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1)
X, Y = np.meshgrid(x, y) 
u = np.zeros((qmax, imax+1, jmax+1))  # 次元を増やすため、jmax+1 の部分を足した。


u[0] = 1 + A * np.sin(2*np.pi*X)
u[1] = A * np.sin(2*np.pi*X)
u[2] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*X))
u[3] = u[2].copy()

fig, ax = plt.subplots(figsize=(6,4), dpi=100)
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0,0,:], label="t = 0")
ims = []
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$y$")
ax2.set_xlabel(r"$x$")
ax2.set_ylabel(r"$u[0]$")

for n in range(nmax):
    u = lax_w_2step_2d(u, calw_fluid_q4v2) # 積分1ステップ進める
    if (n+1) % 100 == 0:
        ax2.plot(x, u[0, 0, :], label=f"t = {dt * (n+1):.2f}")
    txt = ax.text(0, 1.01, f"t = {dt * (n + 1):.1f}", fontsize=15)
    im = ax.contourf(X, Y, u[0, :, :], 20, cmap=plt.cm.viridis)
    ims.append(im.collections+[txt]) # im.collectionsというところに等高線プロットの情報が入っているらしい。これと文字列を結合して、経過時間を表示している。
cbar = fig.colorbar(im)
cbar.set_label(r"Density $\rho$", fontsize=15)
fig.tight_layout()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
# HTML(anim.to_html5_video()) 
ax2.legend()
plt.close(fig)  

calwを変える

$$ \begin{aligned} \begin{gathered} \mathbf{u} \equiv\left(\begin{array}{c} \rho \\ \rho v_{x} \\ \rho v_{y} \\ e \end{array}\right), \mathbf{w}_{1} \equiv\left(\begin{array}{c} \rho v_{x} \\ \rho v_{x}^{2}+p \\ \rho v_{x} v_{y} \\ (e+p) v_{x} \end{array}\right), \mathbf{w}_{2} \equiv\left(\begin{array}{c} \rho v_{y} \\ \rho v_{x} v_{y} \\ \rho v_{y}^{2}+p \\ (e+p) v_{y} \end{array}\right) \end{gathered} \end{aligned} $$

$$ e=\frac{\rho}{2} v^{2}+\frac{p}{\gamma-1} $$

はエネルギー密度

def calw_fluid_2d(u):
    w1 = w2 = np.empty_like(u)
    rho = u[0].copy()
    rho_vx = u[1].copy()
    rho_vy = u[2].copy()
    e = u[3].copy()
    
    vx = rho_vx / rho
    vy = rho_vy / rho
    v2 = vx ** 2 + vy ** 2 # velocity squared
    p = (gamma - 1) * (e - rho * v2 / 2)
    w1[0] = rho_vx
    w1[1] = rho_vx * vx + p
    w1[2] = rho_vx * vy
    w1[3] = (e + p) * vx
    
    w2[0] = rho_vy
    w2[1] = rho_vy * vx
    w2[2] = rho_vy * vy + p
    w2[3] = (e + p) * vy
    
    return w1, w2

imax = 50
jmax = 50
nmax = 200
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 1e-2
kap1 = dt/dx
kap2 = dt/dy
# kap = dt/dx
gamma = 5/3
A = 0.1
fname = "data/courseB-6-A01.dat"
j_flag = 0

x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1)
u = np.zeros((qmax, jmax+1, imax+1))  # u[q,y,x]  次元を増やすため、jmax+1 の部分を足した。
X, Y = np.meshgrid(x, y)

u[0] = 1 + A * np.sin(2*np.pi*X)
u[1] = A * np.sin(2*np.pi*X)
u[2] = u[1].copy()
u[3] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*X))

fig, ax = plt.subplots(figsize=(6,4))
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0, 0,:], label="t = 0")
ims = []
ax.set_xlabel("$x$")
ax.set_ylabel("u[0]")
ax2.set_xlabel("$x$")
ax2.set_ylabel("u[0]")

# テキストファイルが存在していたら削除する
if os.path.exists(fname):
    os.remove(fname)
f = open(fname, "a")  # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
    np.savetxt(f, u[0, 0, :].reshape(1, len(u[1, 0, :])))
    u = lax_w_2step_2d(u, calw_fluid_2d) # 積分1ステップ進める
    if (n+1) % 100 == 0:
        ax2.plot(x, u[0, 0, :], label=f"t = {dt * (n+1):.2f}")
    # animation
    txt = ax.text(0, 1, f"t = {dt * (n + 1):.1f}")
    im, = ax.plot(x, u[0, 0, :], "black")
    ims.append([im, txt])
f.close()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
ax2.legend()
plt.close(fig)  

imax = 50
jmax = 50
nmax = 200
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 1e-2
kap1 = dt/dx
kap2 = dt/dy
gamma = 5/3
A = 0.1

x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1)
X, Y = np.meshgrid(x, y) 
u = np.zeros((qmax, imax+1, jmax+1))  # 次元を増やすため、jmax+1 の部分を足した。

u[0] = 1 + A * np.sin(2*np.pi*X)
u[1] = A * np.sin(2*np.pi*X)
u[2] = u[1].copy()
u[3] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*X))

fig, ax = plt.subplots(figsize=(6,4), dpi=100)
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0, 0 ,:], label="t = 0")
ims = []
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$y$")
ax2.set_xlabel(r"$y$")
ax2.set_ylabel(r"$u[0]$")

for n in range(nmax):
    u = lax_w_2step_2d(u, calw_fluid_2d) # 積分1ステップ進める
    if (n+1) % 100 == 0:
        ax2.plot(x, u[0, 0 ,:], label=f"t = {dt * (n+1):.2f}")
    txt = ax.text(0, 1.01, f"t = {dt * (n + 1):.1f}", fontsize=15)
    im = ax.contourf(X, Y, u[0, : ,:], 20, cmap=plt.cm.viridis)
    ims.append(im.collections+[txt]) # im.collectionsというところに等高線プロットの情報が入っているらしい。これと文字列を結合して、経過時間を表示している。
cbar = fig.colorbar(im)
cbar.set_label(r"Density $\rho$", fontsize=15)
fig.tight_layout()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
# HTML(anim.to_html5_video()) 
ax2.legend()
plt.close(fig)  

 !diff -s data/courseB-6-A01.dat data/courseB-2-A01.dat 

課程5のものとプロットを比較してみよう。 概ね等しいが、僅かに異なっていることがわかる。

f1 = np.loadtxt("data/course5-ex1-A01.dat")
f2 = np.loadtxt("data/courseB-6-A01.dat")
fig, ax = plt.subplots()
ax.plot(x, f1[0,:] - f2[0,:], label="n = 0")
# ax.plot(x, f1[0,:], "--", label="課程5")
# ax.plot(x, f2[0,:], ".", label="付録B")
for i in range(99, 202, 100):
    # ax.plot(x, f1[i,:], "--", label="課程5")
    # ax.plot(x, f2[i,:], ".", label="付録B")
    ax.plot(x, f1[i,:] - f2[i,:], label=f"n = {i}")
ax.legend()

<matplotlib.legend.Legend at 0x13fe164f0>

y方向への伝播を見る。

u[i]が$\text{imax} \times \text{jmax}$行列で、$\mathbf{x} \otimes \mathbf{y}$に等しいので、それの転置を取れば$\mathbf{y} \otimes \mathbf{x}$になり、$x, y$を入れ替えた場合を考えられる。

imax = 50
jmax = 50
nmax = 200
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 1e-2
kap1 = dt/dx
kap2 = dt/dy
gamma = 5/3
A = 0.1

x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1)
X, Y = np.meshgrid(x, y) 
u = np.zeros((qmax, imax+1, jmax+1))  # 次元を増やすため、jmax+1 の部分を足した。

u[0] = 1 + A * np.sin(2*np.pi*Y)
u[1] = A * np.sin(2*np.pi*Y)
u[2] = u[1].copy()
u[3] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*Y))

fig, ax = plt.subplots(figsize=(6,4), dpi=100)
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0,:,0], label="t = 0")
ims = []
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$y$")
ax2.set_xlabel(r"$y$")
ax2.set_ylabel(r"$u[0]$")

for n in range(nmax):
    u = lax_w_2step_2d(u, calw_fluid_2d) # 積分1ステップ進める
    if (n+1) % 100 == 0:
        ax2.plot(x, u[0, :, 0], label=f"t = {dt * (n+1):.2f}")
    txt = ax.text(0, 1.01, f"t = {dt * (n + 1):.1f}", fontsize=15)
    im = ax.contourf(X, Y, u[0, :, :], 20, cmap=plt.cm.viridis)
    ims.append(im.collections+[txt]) # im.collectionsというところに等高線プロットの情報が入っているらしい。これと文字列を結合して、経過時間を表示している。
cbar = fig.colorbar(im)
cbar.set_label(r"Density $\rho$", fontsize=15)
fig.tight_layout()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
# HTML(anim.to_html5_video()) 
ax2.legend()
plt.close(fig)  

animationから明らかなように、y方向への伝播が観測できる。

斜め方向への伝播を見る。

次のコードを実行してわかるように、$t = 0.6$程度で計算が破綻している。 x,y方向への伝播ではみえなかったが、後の爆発のシムレーションの際に左上、右下の領域での計算結果がおかしくなっていることが見えるため、この非対称性によって計算がうまくいっていないのだろう。

計算が対称的でないのはラックスウェンドロフ法の関数に起因するであろうから、見直すべきかもしれない。

imax = 50
jmax = 50
nmax = 70
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 1e-2
kap1 = dt/dx
kap2 = dt/dy
gamma = 5/3
A = 0.1

x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1)
X, Y = np.meshgrid(x, y) 
u = np.zeros((qmax, imax+1, jmax+1))  # 次元を増やすため、jmax+1 の部分を足した。

X_roll = np.zeros_like(X)
for j in range(jmax):
    X_roll[j, :] = np.roll(X[j], -j) 

u[0] = 1 + A * np.sin(2*np.pi*X_roll)
u[1] = A * np.sin(2*np.pi*X_roll)
u[2] = u[1].copy()
u[3] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*X_roll))

# u = np.roll(u, 35, axis=1)

fig, ax = plt.subplots(figsize=(6,4), dpi=100)
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0,0,:], label="t = 0")
ims = []
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$y$")
ax2.set_xlabel(r"$x$")
ax2.set_ylabel(r"$u[0]$")

for n in range(nmax):
    u = lax_w_2step_2d(u, calw_fluid_2d) # 積分1ステップ進める
    if (n+1) % 100 == 0:
        ax2.plot(x, u[0, 0, :], label=f"t = {dt * (n+1):.2f}")
    txt = ax.text(0, 1.01, f"t = {dt * (n + 1):.1f}", fontsize=15)
    im = ax.contourf(X, Y, u[0, :, :], 20, cmap=plt.cm.viridis)
    ims.append(im.collections+[txt]) # im.collectionsというところに等高線プロットの情報が入っているらしい。これと文字列を結合して、経過時間を表示している。
cbar = fig.colorbar(im)
cbar.set_label(r"Density $\rho$", fontsize=15)
fig.tight_layout()
anim = ArtistAnimation(fig, ims, interval=50)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
# HTML(anim.to_html5_video()) 
ax2.legend()
plt.close(fig)  

X_roll = np.zeros_like(X)
for j in range(jmax):
    X_roll[j, :] = np.roll(X[j], j) 
plt.contourf(X, Y, X_roll)
# X[0, :]

<matplotlib.contour.QuadContourSet at 0x1681254c0>

ケルビンヘルムホルツ不安定性

以下、摂動を加える必要あり

レイリーテイラー不安定性 密度を上部、下部で変えて重力を設定。

ケルビンヘルムホルツ不安定性 vxを、下半分は右側、上半分は左側に動くように設定する。

imax = 50
jmax = 50
nmax = 160
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 5e-3
kap1 = dt/dx
kap2 = dt/dy
gamma = 5/3
A = 0.1
nplot = 5

x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1)
X, Y = np.meshgrid(x, y) 
u = np.zeros((qmax, imax+1, jmax+1))  # 次元を増やすため、jmax+1 の部分を足した。

w0 = 0.1
sigma = 0.05/np.sqrt(2.)
gamma = 5/3.
rho = np.ones_like(X)
vx = -0.5 + (np.abs(Y-0.5)<0.25)
vy = w0*np.sin(4*np.pi*X) * ( np.exp(-(Y-0.25)**2/(2 * sigma**2)) + np.exp(-(Y-0.75)**2/(2*sigma**2)) )
vz = 0*X
P = 0*X + 2.5

u[0] = rho
u[1] = rho * vx
u[2] = rho * vy
u[3] = rho * (vx ** 2 + vy ** 2) / 2 + P/(gamma - 1)
# plt.contourf(X, Y, u[1])
# plt.colorbar()

# plt.show()
# plt.contourf(X, Y, u[0])
# plt.colorbar()

# plt.show()

# u[0] = np.tile(1 + A * np.sin(2*np.pi*x), (jmax+1,1)) # density rho
# u[1] = np.tile(A * np.sin(2*np.pi*x), (jmax+1,1))     # momentum rho v 
# u[2] = np.tile(0.9 * (1 + gamma * A * np.sin(2*np.pi*x)), (jmax+1,1)) # energy e
# u[3] = u[2].copy()

# fig, ax = plt.subplots(figsize=(6,4), dpi=100)
fig2, ax2 = plt.subplots(2, 1, figsize=(6,12))
# fig2.colorbar(plot)
# ims = []
# ax.set_xlabel(r"$x$")
# ax.set_ylabel(r"$y$")
# ax2.set_xlabel(r"$y$")
# ax2.set_ylabel(r"$u[0]$")
plot = ax2[0].contourf(X, Y, u[1, :, ])

# txt = ax.text(0, 1.01, f"t = 0", fontsize=15)
# im = ax.contourf(X, Y, u[0, :, :], 20, cmap=plt.cm.viridis)
# ims.append(im.collections+[txt]) # im.collectionsというところに等高線プロットの情報が入っているらしい。これと文字列を結合して、経過時間を表示している。

for n in range(nmax):
    u = lax_w_2step_2d(u, calw_fluid_2d) # 積分1ステップ進める
    # if n == 0:
    #     ax2.contourf(X, Y, u[1, :, :])
    # txt = ax.text(0, 1.01, f"t = {dt * (n + 1):.1f}", fontsize=15)
    # im = ax.contourf(X, Y, u[0, :, :], 20, cmap=plt.cm.viridis)
    # ims.append(im.collections+[txt]) # im.collectionsというところに等高線プロットの情報が入っているらしい。これと文字列を結合して、経過時間を表示している。
plot = ax2[1].contourf(X, Y, u[1, :, ])

# cbar = fig.colorbar(im)
# cbar.set_label(r"Density $\rho$", fontsize=15)
# fig.tight_layout()
# anim = ArtistAnimation(fig, ims, interval=30)
# video = anim.to_html5_video()
# html = display.HTML(video)
# display.display(html)
# ax2.legend()
# plt.close(fig)  

imax = 50
jmax = 50
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 2e-3
t_end = 2
nmax = int(t_end / dt)
kap1 = dt/dx
kap2 = dt/dy
gamma = 5/3
A = 0.1
plot_index = 1
nplot = 10
mabiki = 1

x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1)
X, Y = np.meshgrid(x, y) 
u = np.zeros((qmax, imax+1, jmax+1))  # 次元を増やすため、jmax+1 の部分を足した。

def lax_w_2step_2d_ver2(u, calw):# 境界条件を変える
    # Firts step
    w1, w2 = calw(u)
    u_first = np.ones_like(u)
    u_first[:,:-1, :-1] = u[:,:-1, :-1] - kap1 * (w1[:, 1:, :-1] - w1[:,:-1, :-1]) - kap2 * (w2[:, :-1, 1:] - w2[:,:-1, :-1])
    u_first[:, :, -1] = u_first[:, :, 0].copy() # BdC
    u_first[:, -1, :] = u_first[:, 0, :].copy() # BdC
    # Second step
    w1, w2 = calw(u_first)
    u_second = np.ones_like(u)
    u_second[:, 1:, 1:] = 0.5 * (u[:, 1:, 1:] + u_first[:, 1:, 1:] - kap1 * (w1[:, 1:, 1:] - w1[:, :-1, 1:]) - kap2 * (w2[:, 1:, 1:] - w2[:, 1:, :-1]))
    u_second[:, :, 0] = u_second[:, :, -1].copy() # BdC
    u_second[:, 0, :] = u_second[:, -1, :].copy() # BdC
    return u_second

w0 = 0.1
sigma = 0.05/np.sqrt(2.)
gamma = 5/3
delta = 0.05
rho = np.ones_like(X)
vx = np.tanh(- (np.abs(Y - 0.5) - 0.25) / delta) 
vy = w0*np.sin(2*np.pi*X) * vx * np.tanh(Y / delta)
P = 0*X + 2.5

# w0 = 0.1
# sigma = 0.05/np.sqrt(2.)
# rho = 1. + (np.abs(Y-0.5) < 0.25)
# vx = -0.5 + (np.abs(Y-0.5)<0.25)
# vy = w0*np.sin(4*np.pi*X) * ( np.exp(-(Y-0.25)**2/(2 * sigma**2)) + np.exp(-(Y-0.75)**2/(2*sigma**2)) )
# P = 2.5 * np.ones(X.shape)

u[0] = rho
u[1] = rho * vx
u[2] = rho * vy
u[3] = rho * (vx ** 2 + vy ** 2) / 2 + P/(gamma - 1)

for ui in u:
    plt.contourf(X, Y, ui)
    plt.colorbar()
    plt.show()

def plot_vector():
    # plt.streamplot(X[::mabiki, ::mabiki], Y[::mabiki, ::mabiki], vx[::mabiki, ::mabiki], vy[::mabiki, ::mabiki], density=1.5)
    plt.quiver(X[::mabiki, ::mabiki], Y[::mabiki, ::mabiki], vx[::mabiki, ::mabiki], vy[::mabiki, ::mabiki])
    plt.title(f"n = {n}, t = {n * dt}")
    plt.show()

for n in range(nmax):
    if n % (nmax / nplot) == 0:
        plot_vector()
    u = lax_w_2step_2d_ver2(u, calw_fluid_2d) # 積分1ステップ進める
    vx = u[1] / u[0]
    vy = u[2] / u[0]
plot_vector()

vx = np.tanh(- (np.abs(Y - 0.5) - 0.01) / delta*10)
plt.contourf(X, Y, vx)
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x13bd32af0>

二次元の爆発

imax = 50
jmax = 50
nmax = 200
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 1e-3
kap1 = dt/dx
kap2 = dt/dy
gamma = 5/3
A = 0.1
nindex = 3
nplot = 10

x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1)
X, Y = np.meshgrid(x, y) 
u = np.zeros((qmax, imax+1, jmax+1))  # 次元を増やすため、jmax+1 の部分を足した。

R = np.sqrt((X - 0.5) ** 2 + (Y - 0.5) ** 2)

u[0] = np.ones_like(X)
u[1] = np.zeros_like(X)
u[2] = np.zeros_like(X)
u[3] = 0.9 * np.ones_like(X) + (1 - 2 * R) * (R < 0.5)

for ui in u:
    plt.contourf(X, Y, ui)
    plt.colorbar()
    plt.show()

# fig, ax = plt.subplots(figsize=(5,4), dpi=100)
# fig2, ax2 = plt.subplots(figsize=(6,4))
# ax2.plot(x, u[nindex,0,:], label="t = 0")
# ims = []
# ax.set_xlabel(r"$x$")
# ax.set_ylabel(r"$y$")
# ax.set_aspect('equal', adjustable='box')
# ax2.set_xlabel(r"$x$")
# ax2.set_ylabel(r"$u[0]$")

for n in range(nmax):
    u = lax_w_2step_2d(u, calw_fluid_2d)
    if (n + 1) % (nmax / nplot) == 0:
        fig, ax = plt.subplots(figsize=(5,4))
        ax.set_aspect('equal', adjustable='box')
        ax.contourf(X, Y, u[3])
        fig.suptitle(f"t = {dt * (n+1):.3f}")
        fig.tight_layout()
        fig.show()

# for n in range(nmax):
#     u = lax_w_2step_2d(u, calw_fluid_2d) # 積分1ステップ進める
#     # if (n+1) % 100 == 0:
#     #     ax2.plot(x, u[nindex, 0, :], label=f"t = {dt * (n+1):.2f}")
#     txt = ax.text(0, 1.01, f"t = {dt * (n + 1):.1f}", fontsize=15)
#     im = ax.contourf(X, Y, u[nindex, :, :], 20, cmap=plt.cm.viridis)
#     ims.append(im.collections+[txt]) # im.collectionsというところに等高線プロットの情報が入っているらしい。これと文字列を結合して、経過時間を表示している。
# cbar = fig.colorbar(im)
# cbar.set_label(r"Energy", fontsize=15)
# fig.tight_layout()
# anim = ArtistAnimation(fig, ims, interval=100)
# video = anim.to_html5_video()
# html = display.HTML(video)
# display.display(html)
# # HTML(anim.to_html5_video()) 
# # ax2.legend()
# plt.close(fig)  

/var/folders/y5/1mt0b_g15zn1vv8_73m7hlb80000gn/T/ipykernel_79809/2864887905.py:50: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.
  fig.show()

u[3] = 0.9 * np.ones_like(X) + (1 - 2 * R) * (R < 0.5)
plt.contourf(X, Y, u[3])
plt.colorbar()
plt.show()

plt.plot(x, u[3,:, int(jmax-10)])

[<matplotlib.lines.Line2D at 0x13beb0bb0>]

付録C プラズマの数値シミュレーション

$(x, y) -> (r, \theta)$

$$ x = r\cos \theta \\ y = r\sin \theta $$$$ e = 0.9 (r > 0.5)\\ e = 0.9 + 1 - 2*x (r < 0.5) $$

方程式

$$ \frac{\partial \mathbf{u}}{\partial t}=-\frac{\partial \mathbf{w}}{\partial x}\\ \mathbf{u} \equiv\left(\begin{array}{c}\rho \\ \rho v_{x} \\ \rho v_{y} \\ \rho v_{z} \\ e \\ B_{y} \\ B_{z}\end{array}\right), \quad \mathbf{w} \equiv\left(\begin{array}{c}\rho v_{x} \\ \rho v_{x}^{2}+p+B^{2} / 2 \\ \rho v_{x} v_{y}-B_{x} B_{y} \\ \rho v_{x} v_{z}-B_{x} B_{z} \\ \left(e+p+B^{2} / 2\right) v_{x}-\left(v_{x} B_{x}+v_{y} B_{y}+v_{z} B_{z}\right) B_{x} \\ -B_{x} v_{y}+v_{x} B_{y} \\ -B_{x} v_{z}+v_{x} B_{z}\end{array}\right) $$

$$ e=\frac{\rho}{2} v^{2} + \frac{B^2}{2} +\frac{p}{\gamma-1} $$

はエネルギー密度